

% Related parameters (will need in natural wastage region)
ssaux = saux - (1/(1-alpha));
 %
 % Make grid of values of unknown constant of integration:
Ahi  = .99*ssaux*(mL^-ssaux)/((mM/mL)^ssaux -1);
Alo  = .01*ssaux*(mL^-ssaux)/((mM/mL)^ssaux -1);
obs_const= 100; 
constvec = linspace(Alo, Ahi, obs_const)';
 %
 %
%% Natural wastage
Gf_mM = (constvec / ssaux) *(mL^ssaux) *((mM/mL)^ssaux -1);	
 %
 %
%% Full replacement 
const_fr = (mM^saux) * constvec;    % This is saux, not ssaux.
Gf_mR    = Gf_mM + (const_fr/(1/(1-alpha)))*(mM^(-1/(1-alpha)) - mR^(-1/(1-alpha)));
gf_mR    = const_fr*(mR^(-1-(1/(1-alpha))));
 %
 %
%% Net Expansion
%
% We will need boundary value of quit rate:
delta_mR = s*lambda  / (1 + aux*s*lambda*log(mR/mM));
auxconst = ( ((1-omega1)/(alpha*c))*mR - (r+delta_mR + ((omega0+r*C)/c)) )*(mR^(-1/(1-alpha)));
%
% density:
gf_ne = zeros(obs_ne, obs_const);
for i = 1:obs_ne
    [grid2_ne, wghts2_ne] = legendre(mR, grid_ne(i), obs_ne);
    %
    % We will need the quit rate in this region:
    delta_ne = ((1-omega1)/(alpha*c))*grid2_ne - (r+((omega0+r*C)/c)) ...
             - auxconst*(grid2_ne.^(1/(1-alpha)));
    %
    % We will need the hiring rate, too.
    % eta(m) = - delta'(m) * q(m) / q'(m);
    % q(m) / q'(m) = (\psi + (1-\psi)*G_worker(m)) / (1-\psi)*g_worker(m)
    %
    Gw_ne = interp1(grid_ne, G_ne, grid2_ne,'linear','extrap');     % Distribution of m over workers.
    gw_ne = interp1(grid_ne, g_ne, grid2_ne,'linear','extrap');
    %
    deltaprime_ne = ((1-omega1)/(alpha*c)) - auxconst*(1/(1-alpha))*(grid2_ne.^(alpha/(1-alpha)));
    qtoqprime_ne  = (ushr + (1-ushr)*Gw_ne) ./ ((1-ushr)*gw_ne);
    eta_ne        = -deltaprime_ne .* qtoqprime_ne;
    nu            = mu - (1-alpha)*(eta_ne - delta_ne);
    %
    %
    rhsint = wghts2_ne'*( (nu-sigma^2)./grid2_ne );
    rhs    = exp( (2/(sigma^2)) * rhsint );
    gf_ne(i,:) = rhs*gf_mR';
end
%
% Solve for unknown constant:
Gf_mU  = wghts_ne'*gf_ne;       % 1 by obs_const
integral_of_g = Gf_mR + Gf_mU';
if min(integral_of_g) >1 || max(integral_of_g) < 1
    disp('ERROR: density cannot integrate to 1');
    pause;
end
% 
% Solution for unknown constant
constsol = interp1(integral_of_g, constvec, 1);
%
% Final density
gf_nw = constsol*(grid_nw.^(ssaux-1)); 
gf_fr = (mM^saux)*constsol*(grid_fr.^(-1-(1/(1-alpha))));
gf_ne = interp1(constvec, gf_ne', constsol);	
gf_ne = gf_ne';	
gf_all= [gf_nw; gf_fr; gf_ne];
if min(gf_ne) <1e-9
    disp('ERROR: Density fnc. in NE region is negative');
    pause;
end
%
% Final c.d.f.
Gf_mM  = interp1(constvec,Gf_mM,constsol);
Gf_mR  = interp1(constvec,Gf_mR,constsol);
Gf_nw  = cumsum( wghts_nw.*gf_nw );               
Gf_fr  = Gf_mM + cumsum( wghts_fr.*gf_fr);         
Gf_ne  = Gf_mR + cumsum( wghts_ne.*gf_ne);    
Gf_all = [Gf_nw; Gf_fr; Gf_ne];
gridf_all = grid_all;
